﻿using System;
using System.Windows.Forms;
using System.IO;

namespace RSProject
{
    public class funcStatistics
    {
        // 平均值
        public void Mean(string P)
        {
            string s_data = P;
            string s_metadata = P + ".hdr";
            int[,] tempArray = new int[1000, 1000];
            byte[,] tempData = new byte[1000, 1000];

            //初始化
            for (int j = 0; j < 1000; j++)
                for (int k = 0; k < 1000; k++)
                    tempArray[j, k] = 0;

            //计算均值
            for (int j = 1; j < Global.head.lines - 1; j++)
            {
                for (int k = 1; k < Global.head.samples - 1; k++)
                {
                    for (int i = 0; i < Global.head.bands; i++)
                    {
                        tempArray[j, k] += Global.intArray[i, j, k];
                    }
                    tempArray[j, k] /= Global.head.bands;
                }
            }

            //写入文件
            Head tempHead = Global.head;
            tempHead.bands = 1;
            Global.data.write_metadata(s_metadata, "bsq",tempHead);
            BinaryWriter bw = new BinaryWriter(File.Open(s_data, FileMode.Create));

            for (int j = 0; j < 1000; j++)
                for (int k = 0; k < 1000; k++)
                    tempData[j, k] = 0;

            for (int j = 0; j < Global.head.lines; j++)
            {
                for (int k = 0; k < Global.head.samples; k++)
                {
                    tempData[j, k] = Convert.ToByte(tempArray[j, k]);
                    bw.Write(tempData[j, k]);
                }
            }
            bw.Close();
            MessageBox.Show("成功！");
        }

        // 标准差
        public void StandardDeviation(string P)
        {
            string s_data = P;
            string s_metadata = P + ".hdr";
            int[,] tempArray = new int[1000, 1000];
            byte[,] tempData = new byte[1000, 1000];

            //初始化
            for (int j = 0; j < 1000; j++)
                for (int k = 0; k < 1000; k++)
                    tempArray[j, k] = 0;

            //计算标准差
            for (int j = 1; j < Global.head.lines - 1; j++)
            {
                for (int k = 1; k < Global.head.samples - 1; k++)
                {
                    double ex = 0, ex2 = 0;
                    for (int i = 0; i < Global.head.bands; i++)
                    {
                        ex += Global.intArray[i, j, k];
                        ex2 += Math.Pow(Global.intArray[i, j, k], 2);
                    }
                    ex /= 9;
                    ex2 /= 9;

                    tempArray[j, k] = (int)Math.Sqrt(ex2 - Math.Pow(ex, 2));
                }
            }

            //写入文件
            Head tempHead = Global.head;
            tempHead.bands = 1;
            Global.data.write_metadata(s_metadata, "bsq", tempHead);
            BinaryWriter bw = new BinaryWriter(File.Open(s_data, FileMode.Create));

            for (int j = 0; j < 1000; j++)
                for (int k = 0; k < 1000; k++)
                    tempData[j, k] = 0;

            for (int j = 0; j < Global.head.lines; j++)
            {
                for (int k = 0; k < Global.head.samples; k++)
                {
                    tempData[j, k] = Convert.ToByte(tempArray[j, k]);
                    bw.Write(tempData[j, k]);
                }
            }
            bw.Close();
            MessageBox.Show("成功！");
        }

        // 单波段均值
        public double[] MeanSingle()
        {
            double[] ex = new double[Global.head.bands];

            for (int i = 0; i < Global.head.bands; i++)
            {
                ex[i] = 0;
            }

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        ex[i] += Global.intArray[i, j, k];
                    }
                ex[i] /= (Global.head.lines * Global.head.samples);
            }

            return ex;
        }

        // 单波段标准差
        public double[] StandardSingle()
        {
            double[] ex = new double[Global.head.bands];
            double[] ex2 = new double[Global.head.bands];
            double[] dx = new double[Global.head.bands];

            for (int i = 0; i < Global.head.bands; i++)
            {
                ex[i] = 0;
                ex2[i] = 0;
                dx[i] = 0;
            }

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        ex[i] += Global.intArray[i, j, k];
                        ex2[i] += Math.Pow(Global.intArray[i, j, k], 2);
                    }
                ex[i] /= (Global.head.lines * Global.head.samples);
                ex2[i] /= (Global.head.lines * Global.head.samples);
                dx[i] = Math.Sqrt(ex2[i] - Math.Pow(ex[i], 2));
            }

            return dx;
        }

        // 协方差矩阵
        public double[,] Covariance()
        {
            double[,] Cov = new double[Global.head.bands, Global.head.bands];
            double[] ex = new double[Global.head.bands];

            //初始化
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.bands; j++)
                {
                    Cov[i, j] = 0;
                }
                ex[i] = 0;
            }

            //计算协方差
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        ex[i] += Global.intArray[i, j, k];
                    }
                ex[i] /= (Global.head.lines * Global.head.samples);
            }

            for (int m = 0; m < Global.head.bands; m++)
            {
                for (int n = 0; n < Global.head.bands; n++)
                {
                    for (int j = 0; j < Global.head.lines; j++)
                    {
                        for (int k = 0; k < Global.head.samples; k++)
                        {
                            Cov[m, n] += (Global.intArray[m, j, k] - ex[m]) * (Global.intArray[n, j, k] - ex[n]);
                        }
                    }
                    Cov[m, n] /= (Global.head.lines * Global.head.samples);
                }
            }

            //显示协方差矩阵
            return Cov;
        }

        // 相关系数矩阵
        public double[,] CorrelationCoefficient()
        {
            double[,] Cov = new double[Global.head.bands, Global.head.bands];
            double[] ex = new double[Global.head.bands];
            double[] ex2 = new double[Global.head.bands];
            double[] dx = new double[Global.head.bands];//这里指标准差
            double[,] Cor = new double[Global.head.bands, Global.head.bands];//相关系数

            //初始化
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.bands; j++)
                {
                    Cov[i, j] = 0;
                }
                ex[i] = 0;
                ex2[i] = 0;
            }

            //计算标准差
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        ex[i] += Global.intArray[i, j, k];
                        ex2[i] += Math.Pow(Global.intArray[i, j, k], 2);
                    }
                ex[i] /= (Global.head.lines * Global.head.samples);
                ex2[i] /= (Global.head.lines * Global.head.samples);
                dx[i] = Math.Sqrt(ex2[i] - Math.Pow(ex[i], 2));
            }

            //计算协方差
            for (int m = 0; m < Global.head.bands; m++)
            {
                for (int n = 0; n < Global.head.bands; n++)
                {
                    for (int j = 0; j < Global.head.lines; j++)
                    {
                        for (int k = 0; k < Global.head.samples; k++)
                        {
                            Cov[m, n] += (Global.intArray[m, j, k] - ex[m]) * (Global.intArray[n, j, k] - ex[n]);
                        }
                    }
                    Cov[m, n] /= (Global.head.lines * Global.head.samples);
                }
            }

            //计算相关系数
            for (int i = 0; i < Global.head.bands; i++)
                for (int j = 0; j < Global.head.bands; j++)
                {
                    Cor[i, j] = Cov[i, j] / (dx[i] * dx[j]);
                }

            return Cor;
        }

        // 极差纹理
        public void Extreme(string P)
        {
            string s_data = P;
            string s_metadata = P + ".hdr";

            int[,,] tempArray = new int[10, 1000, 1000];
            byte[,,] dataArray = new byte[100, 1000, 1000];

            for (int i = 0; i < 10; i++)
                for (int j = 0; j < 1000; j++)
                    for (int k = 0; k < 1000; k++)
                        tempArray[i, j, k] = 0;

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 1; j < Global.head.lines - 1; j++)
                {
                    for (int k = 1; k < Global.head.samples - 1; k++)
                    {
                        int max = Global.intArray[i, j, k], min = Global.intArray[i, j, k];
                        for (int m = j - 1; m <= j + 1; m++)
                        {
                            for (int n = k - 1; n <= k + 1; n++)
                            {
                                if (Global.intArray[i, m, n] > max)
                                    max = Global.intArray[i, m, n];
                                if (Global.intArray[i, m, n] < min)
                                    min = Global.intArray[i, m, n];
                            }
                        }
                        tempArray[i, j, k] = max - min;
                    }
                }
            }

            //写入文件
            Global.data.write_metadata(s_metadata, "bsq",Global.head);
            BinaryWriter bw = new BinaryWriter(File.Open(s_data, FileMode.Create));

            for (int i = 0; i < 100; i++)
                for (int j = 0; j < 1000; j++)
                    for (int k = 0; k < 1000; k++)
                        dataArray[i, j, k] = 0;

            for (int i = 0; i < Global.head.bands; i++)
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        dataArray[i, j, k] = Convert.ToByte(tempArray[i, j, k]);
                    }

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                {
                    for (int k = 0; k < Global.head.samples; k++)
                        bw.Write(dataArray[i, j, k]);
                }
            }
            bw.Close();
            MessageBox.Show("成功！");
        }

        // 计算直方图
        public void funcHistogram()
        {
            int band;
            Bands hist = new Bands(Global.head.bands, 1);

            // 选择波段
            if (hist.ShowDialog() == DialogResult.OK)
            {
                band = hist.band1;
                dataView view = new dataView(band, -1);
                view.ShowDialog();
            }
        }

        // 计算累计直方图
        public void CumuHistogram()
        {
            int band;
            Bands hist = new Bands(Global.head.bands, 1);

            // 选择波段
            if (hist.ShowDialog() == DialogResult.OK)
            {
                band = hist.band1;
                dataView view = new dataView(band, -2);
                view.ShowDialog();
            }
        }

        // 计算联合直方图
        public void JointHistogram(ref int band1, ref int band2)
        {
            Bands hist = new Bands(Global.head.bands, 2);

            // 选择波段
            if (hist.ShowDialog() == DialogResult.OK)
            {
                band1 = hist.band1;
                band2 = hist.band2;
            }
        }

        // 计算共生矩阵
        public void SymbioticMatrix()
        {
            int a, b;
            SymbioticMatrix sm = new SymbioticMatrix();
            if (sm.ShowDialog() == DialogResult.OK)
            {
                Global.F = new int[256, 256];
                a = sm.a;
                b = sm.b;
                for (int u = 0; u <= 255; u++)
                {
                    for (int v = 0; v <= 255; v++)
                    {
                        Global.F[u, v] = 0;
                    }
                }

                int u0, v0;

                for (int i = 0; i < Global.head.bands; i++)
                {
                    for (int j = 0; j < Global.head.lines; j++)
                    {
                        for (int k = 0; k < Global.head.samples; k++)
                        {
                            u0 = Global.intArray[i, j, k];
                            v0 = Global.intArray[i, j + b, k + a];
                            Global.F[u0, v0]++;
                        }
                    }
                }
            }
        }

        // 计算对比度
        public int Contrast()
        {
            int contrast = 0;
            SymbioticMatrix();
            for (int u = 0; u <= 255; u++)
            {
                for (int v = 0; v <= 255; v++)
                {
                    contrast += (u - v) * (u - v) * Global.F[u, v];
                }
            }
            return contrast;
        }

        // 计算频率
        public int[,] Frequency(int[,,] intdata,Head head)
        {
            int[,] frestr = new int[10, 256];

            for (int i = 0; i < head.bands; i++)
            {
                for (int j = 0; j < head.lines; j++)
                {
                    for (int k = 0; k < head.samples; k++)
                    {
                        frestr[i, intdata[i, j, k]]++;
                    }
                }
            }
            return frestr;
        }
    }
}
